#This is to get:
#Table1
#TableA3
#TableA5
#TableA8
#TableA9
#TableA10

rm(list = ls(all.names = TRUE))
library("stargazer")
library("haven")
library("lfe")
library("rgdal")
library("margins")
library("broom")
library("spdep")
library("multiwayvcov")
library("lmtest")
library("rgeos")
library("spatialreg")

setwd("/Users/bgpopescu/Dropbox/Legacies_Central_Europe/")
setwd("C:/Users/bogdanp/Dropbox/Legacies_Central_Europe/")


#Step0: Function for calculating Robust SE
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, mynewspdf2@data$ctrnum, force_posdef=T))[, 2]
  return(res)}

#Step1: Reading DTA files
nuts_df<-read_dta("./data/nuts_data.dta")

#Step2: Reading shapefile for NUTS polygons (counties)
nuts_sh <- readOGR(dsn="./data/shape_files.gdb",
                layer="nuts_NONEU_and_EU", use_iconv = TRUE, encoding = "UTF-8")

#Step3: Merge
nuts<-merge(nuts_sh, nuts_df, by.x = "MERGE_ID", by.y = "merge_id")
nuts2<-subset(nuts, otthist==1)
rm(nuts_df, nuts_sh, nuts)

#Step4: Coding communist countries
nuts2@data$communist<-ifelse(nuts2@data$ctr=='ES' | nuts2@data$ctr=='IT' 
                             | nuts2@data$ctr=='EL' | nuts2@data$ctr=='CY', 0, 1)
unique(nuts2@data$communist)

#########
#MODEL 1#
#########

#Step_M1_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$p_ottoman) &
                  !is.na(nuts2@data$lngdpcap)),]

#Step_M1_2: Create list of neighbors
#The function "poly2nb" builds a neighbours list based on 
#regions with contiguous boundaries, that is sharing one or more boundary points.
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M1_3: Convert nb to listw type
#The "nb2listw" function supplements a neighbours list with spatial weights for the chosen coding scheme. 
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M1_4: Run regression
felm1_res<-felm(lngdpcap~p_ottoman
               |0|0|ctrnum,
             data=mynewspdf2@data)
summary(felm1_res)
#Saving Cluster SE
felm1_res_cse<-felm1_res$cse

#Step_M1_5: Run regression with distance to Istanbul
felm1_res_di<-felm(lngdpcap~p_ottoman + istanbul_distance
                |0|0|ctrnum,
                data=mynewspdf2@data)
summary(felm1_res_di)
#Saving Cluster SE
felm1_res_cse_di<-felm1_res_di$cse


#Step_M1_6: Run regression with distance to Vienna
felm1_res_vi<-felm(lngdpcap~p_ottoman + vienna_distance
                   |0|0|ctrnum,
                   data=mynewspdf2@data)
summary(felm1_res_vi)
#Saving Cluster SE
felm1_res_cse_vi<-felm1_res_vi$cse

#Step_M1_7: Spatial Filtering
#SpatialFiltering selects eigenvectors in a semi-parametric spatial filtering approach to removing spatial 
#dependence from linear models
#generates Markov Chain Monte Carlo (MCM) and calculates its eigenvectors and 
#then uses a simultaneous autoregressive model to judiciously select a subset of the spatial patterns
#Spatial filtering chooses eigenvectors that reduce the residual spatial autocorrelation in 
#the error of the model with covariates. 
#The lag form adds the covariates in assessment of which eigenvectors to choose, 
#but does not use them in constructing the eigenvectors.

#Note that ExactEV is False. Set ExactEV=TRUE to use exact expectations and variances rather 
#than the expectation and variance of Moran's I from the previous iteration, default FALSE
#When ExactEV is TRUE, it takes 25 mins
lm1_res_sp<-SpatialFiltering(lngdpcap~p_ottoman,
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M1_8: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm1_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M1_9: Include fitted results
x<-as.formula(paste("lngdpcap~p_ottoman + ", paste(eigen_string), "|0|0|ctrnum"))
felm1_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm1_res_fitted)
#Saving Cluster SE
felm1_res_fitted_cse<-felm1_res_fitted$cse

#Step_M1_10: Save Moran I stats
lm1_moran<-moran.test(x=resid(felm1_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm1_moran_p<-moran.test(x=resid(felm1_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm1_moran_star<-paste(round(lm1_moran, digits = 3), "",
                        ifelse(lm1_moran_p<.01,"***",
                               ifelse(lm1_moran_p<.05,"**",
                                      ifelse(lm1_moran_p<.1,"*",
                                             ""))), "", "", sep = "")


#########
#MODEL 2#
#########

#Step_M2_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$p_ottoman) &
                              !is.na(nuts2@data$lngdpcap)),]

#Step_M2_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$MERGE_ID)

#Step_M2_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M2_4: Run regression
felm2_res<-felm(lngdpcap~p_ottoman 
                |ctrnum|0|ctrnum,
                data=mynewspdf2@data)
summary(felm2_res)
#Saving Cluster SE
felm2_res_cse<-felm2_res$cse

#Step_M2_5: Run regression with distance to Istanbul
felm2_res_di<-felm(lngdpcap~p_ottoman + istanbul_distance 
                   |ctrnum |0|ctrnum,
               data=mynewspdf2@data)
summary(felm2_res_di)
#Saving Cluster SE
felm2_res_cse_di<-felm2_res_di$cse

#Step_M2_6: Run regression with distance to Vienna
felm2_res_vi<-felm(lngdpcap~p_ottoman + vienna_distance 
                   |ctrnum |0|ctrnum,
                   data=mynewspdf2@data)
summary(felm2_res_vi)
#Saving Cluster SE
felm2_res_cse_vi<-felm2_res_vi$cse

#Step_M2_7: Spatial Filtering
lm2_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M2_8: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm2_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M2_9: Include fitted results
x<-as.formula(paste("lngdpcap~p_ottoman + ", paste(eigen_string), "|ctrnum|0|ctrnum"))
felm2_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm2_res_fitted)
#Saving Cluster SE
felm2_res_fitted_cse<-felm2_res_fitted$cse

#Step_M2_10: Save Moran I stats
lm2_moran<-moran.test(x=resid(felm2_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm2_moran_p<-moran.test(x=resid(felm2_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm2_moran_star<-paste(round(lm2_moran, digits = 3), "",
                        ifelse(lm2_moran_p<.01,"***",
                               ifelse(lm2_moran_p<.05,"**",
                                      ifelse(lm2_moran_p<.1,"*",
                                             ""))), "", "", sep = "")


#########
#MODEL 3#
#########

#Step_M3_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                            !is.na(nuts2@data$p_ottoman) &
                            !is.na(nuts2@data$ctrnum) &
                            !is.na(nuts2@data$tavg) &
                            !is.na(nuts2@data$lograin) &
                            !is.na(nuts2@data$logelev) &
                            !is.na(nuts2@data$logslope) & 
                            !is.na(nuts2@data$wheat_rainfed) &
                            !is.na(nuts2@data$log_sea_lines_distance) & 
                            !is.na(nuts2@data$log_large_rivers_distance) & 
                            !is.na(nuts2@data$point_x) & 
                            !is.na(nuts2@data$point_y)),]

#Step_M3_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M3_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M3_4: Run regression
felm3_res<-felm(lngdpcap~p_ottoman +
              tavg + lograin + logelev +
              logslope + wheat_rainfed + log_sea_lines_distance +
              log_large_rivers_distance + point_x + point_y
              |ctrnum | 0 | ctrnum,
            data=mynewspdf2@data)
summary(felm3_res)
#Saving Cluster SE
felm3_res_cse<-felm3_res$cse

#Step_M3_5: Run regression with distance to Istanbul
felm3_res_di<-felm(lngdpcap~p_ottoman + tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y+
                 istanbul_distance 
                 |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm3_res_di)
#Saving Cluster SE
felm3_res_cse_di<-felm3_res_di$cse

#Step_M3_6: Run regression with distance to Vienna
felm3_res_vi<-felm(lngdpcap~p_ottoman + tavg + lograin + logelev +
                     logslope + wheat_rainfed + log_sea_lines_distance +
                     log_large_rivers_distance + point_x + point_y+
                     vienna_distance 
                   |ctrnum | 0 | ctrnum,
                   data=mynewspdf2@data)
summary(felm3_res_vi)
#Saving Cluster SE
felm3_res_cse_vi<-felm3_res_vi$cse

#Step_M3_7: Spatial Filtering
lm3_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y + factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, 
                              zero.policy = T)

#Step_M3_8: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm3_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M3_9: Include fitted results
x<-as.formula(paste("lngdpcap~p_ottoman + tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y + ", paste(eigen_string), 
                    "|ctrnum | 0 | ctrnum"))

felm3_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm3_res_fitted)
#Saving Cluster SE
felm3_res_fitted_cse<-felm3_res_fitted$cse

#Step_M3_10: Save Moran I stats
lm3_moran<-moran.test(x=resid(felm3_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm3_moran_p<-moran.test(x=resid(felm3_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm3_moran_star<-paste(round(lm3_moran, digits = 3), "",
                       ifelse(lm3_moran_p<.01,"***",
                              ifelse(lm3_moran_p<.05,"**",
                                     ifelse(lm3_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 4#
#########

#Step_M4_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                              !is.na(nuts2@data$p_ottoman) &
                              !is.na(nuts2@data$route_density_2) &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]

#Step_M4_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M4_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M4_4: Run regression
felm4_res<-felm(lngdpcap~p_ottoman + route_density_2 +
              tavg + lograin + logelev +
              logslope + wheat_rainfed + log_sea_lines_distance +
              log_large_rivers_distance + point_x + point_y 
              |ctrnum | 0 | ctrnum,
            data=mynewspdf2@data)
summary(felm4_res)
#Saving Cluster SE
felm4_res_cse<-felm4_res$cse

#Step_M4_5: Run regression with distance to Istanbul
felm4_res_di<-felm(lngdpcap~p_ottoman + route_density_2 +
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y +
                 istanbul_distance 
                |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm4_res_di)
#Saving Cluster SE
felm4_res_cse_di<-felm4_res_di$cse

#Step_M4_6: Run regression with distance to Vienna
felm4_res_vi<-felm(lngdpcap~p_ottoman + route_density_2 +
                     tavg + lograin + logelev +
                     logslope + wheat_rainfed + log_sea_lines_distance +
                     log_large_rivers_distance + point_x + point_y +
                     vienna_distance 
                   |ctrnum | 0 | ctrnum,
                   data=mynewspdf2@data)
summary(felm4_res_vi)
#Saving Cluster SE
felm4_res_cse_vi<-felm4_res_vi$cse


#Step_M4_7: Spatial Filtering
lm4_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + route_density_2 +
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y+
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M4_8: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm4_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)


#Step_M4_9: Include fitted results
x<-as.formula(paste("lngdpcap~p_ottoman + route_density_2 + tavg + lograin + 
                        logelev + logslope + wheat_rainfed + log_sea_lines_distance +
                        log_large_rivers_distance + point_x + point_y  +", paste(eigen_string),
                       "|ctrnum | 0 | ctrnum"))

felm4_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm4_res_fitted)
#Saving Cluster SE
felm4_res_fitted_cse<-felm4_res_fitted$cse

#Step_M4_10: Save Moran I stats
lm4_moran<-moran.test(x=resid(felm4_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm4_moran_p<-moran.test(x=resid(felm4_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm4_moran_star<-paste(round(lm4_moran, digits = 3), "",
                       ifelse(lm4_moran_p<.01,"***",
                              ifelse(lm4_moran_p<.05,"**",
                                     ifelse(lm4_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 5#
#########

#Step_M5_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                              !is.na(nuts2@data$p_ottoman) &
                              !is.na(nuts2@data$yearpos) &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]


#Step_M5_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M5_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M5_4: Run regression
lm5_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(lm5_res)


#Step_M5_5: Saving coefficients and SE
coef_lm5_res<-summary(margins(lm5_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$AME
se_lm5_res<-summary(margins(lm5_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$SE

d5 <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 2), nc = 2))
names(d5) <- c("lngdpcap", "p_ottoman")
f5 <- as.formula("lngdpcap ~ 0 + p_ottoman")
p5 <- lm(f5, d5)
p5$coefficients<-coef_lm5_res
p5_se<- setNames(se_lm5_res, c("p_ottoman"))
p5$df.residual<-NA


#Step_M5_6: Run regression with distance to Istanbul
lm5_res_di<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum) +
               istanbul_distance,
             data=mynewspdf2@data)
summary(lm5_res_di)

#Step_M5_7: Saving coefficients and SE
coef_lm5_res_di<-summary(margins(lm5_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1650)))$AME
se_lm5_res_di<-summary(margins(lm5_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1650)))$SE

d5_di <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d5_di) <- c("lngdpcap", "p_ottoman", "istanbul_distance")
f5_di <- as.formula("lngdpcap ~ 0 + p_ottoman + istanbul_distance")
p5_di <- lm(f5_di, d5_di)
p5_di$coefficients<-coef_lm5_res_di
p5_se_di<- setNames(se_lm5_res_di, c("p_ottoman", "istanbul_distance"))
p5_di$df.residual<-NA

#Step_M5_8: Run regression with distance to Vienna
lm5_res_vi<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  vienna_distance+
                 factor(ctrnum),
                data=mynewspdf2@data)
summary(lm5_res_vi)

#Step_M5_9: Saving coefficients and SE
coef_lm5_res_vi<-summary(margins(lm5_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1650)))$AME
se_lm5_res_vi<-summary(margins(lm5_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1650)))$SE

d5_vi <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d5_vi) <- c("lngdpcap", "p_ottoman", "vienna_distance")
f5_vi <- as.formula("lngdpcap ~ 0 + p_ottoman + vienna_distance")
p5_vi <- lm(f5_vi, d5_vi)
p5_vi$coefficients<-coef_lm5_res_vi
p5_se_vi<- setNames(se_lm5_res_vi, c("p_ottoman", "vienna_distance"))
p5_vi$df.residual<-NA

#Step_M5_10: Run regression with Trade Route Density
lm5_res_tr<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  route_density_2 +
                 factor(ctrnum),
                data=mynewspdf2@data)
summary(lm5_res_tr)

#Step_M5_11: Saving coefficients and SE
coef_lm5_res_tr<-summary(margins(lm5_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1650)))$AME
se_lm5_res_tr<-summary(margins(lm5_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1650)))$SE

d5_tr <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d5_tr) <- c("lngdpcap", "p_ottoman", "route_density_2")
f5_tr <- as.formula("lngdpcap ~ 0 + p_ottoman + route_density_2")
p5_tr <- lm(f5_tr, d5_tr)
p5_tr$coefficients<-coef_lm5_res_tr
p5_se_tr<- setNames(se_lm5_res_tr, c("p_ottoman", "route_density_2"))
p5_tr$df.residual<-NA

#Step_M5_12: Spatial Filtering
lm5_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + yearpos+ 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y +
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M5_13: Creating the strings of eigenvalues
mynewspdf3<-mynewspdf2@data
mynewspdf4<-cbind(mynewspdf3, fitted(lm5_res_sp))
eigen_vars<-c(names(mynewspdf4[,grepl( "vec" , names(mynewspdf4))]))

#Step_M5_14: Include fitted results
lm5_res_fitted<-lm(as.formula(
               paste("lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum)+", paste(eigen_vars, collapse = "+"))),
             data=mynewspdf4)
summary(lm5_res_fitted)

#Step_M5_15: Saving coefficients and SE
coef_lm5_res_fitted<-summary(margins(lm5_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$AME
se_lm5_res_fitted<-summary(margins(lm5_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1650)))$SE

d5_fitted <- as.data.frame(matrix(rnorm(nrow(mynewspdf4) * 2), nc = 2))
names(d5_fitted) <- c("lngdpcap", "p_ottoman")
f5_fitted <- as.formula("lngdpcap ~ 0 + p_ottoman")
p5_fitted <- lm(f5_fitted, d5_fitted)
p5_fitted$coefficients<-coef_lm5_res_fitted
p5_fitted_se<- setNames(se_lm5_res_fitted, c("p_ottoman"))
p5_fitted$df.residual<-NA

#Step_M5_16: Save Moran I stats
lm5_moran<-moran.test(x=resid(lm5_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm5_moran_p<-moran.test(x=resid(lm5_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm5_moran_star<-paste(round(lm5_moran, digits = 3), "",
                       ifelse(lm5_moran_p<.01,"***",
                              ifelse(lm5_moran_p<.05,"**",
                                     ifelse(lm5_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 6#
#########

#Step_M6_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                              !is.na(nuts2@data$p_ottoman) &
                              !is.na(nuts2@data$yearpos) &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]


#Step_M6_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M6_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M6_4: Run regression
lm6_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(lm6_res)

#Step_M6_5: Saving coefficients and SE
coef_lm6_res<-summary(margins(lm6_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$AME
se_lm6_res<-summary(margins(lm6_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$SE

d6 <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 2), nc = 2))
names(d6) <- c("lngdpcap", "p_ottoman")
f6 <- as.formula("lngdpcap ~ 0 + p_ottoman")
p6 <- lm(f6, d6)
p6$coefficients<-coef_lm6_res
p6_se<- setNames(se_lm6_res, c("p_ottoman"))
p6$df.residual<-NA


#Step_M6_6: Run regression with distance to Istanbul
lm6_res_di<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum) +
               istanbul_distance,
             data=mynewspdf2@data)
summary(lm6_res_di)

#Step_M6_7: Saving coefficients and SE
coef_lm6_res_di<-summary(margins(lm6_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1700)))$AME
se_lm6_res_di<-summary(margins(lm6_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1700)))$SE

d6_di <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d6_di) <- c("lngdpcap", "p_ottoman", "istanbul_distance")
f6_di <- as.formula("lngdpcap~ 0 + p_ottoman + istanbul_distance")
p6_di <- lm(f6_di, d6_di)
p6_di$coefficients<-coef_lm6_res_di
p6_se_di<- setNames(se_lm6_res_di, c("p_ottoman", "istanbul_distance"))
p6_di$df.residual<-NA


#Step_M6_8: Run regression with distance to Vienna
lm6_res_vi<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  factor(ctrnum) +
                  vienna_distance,
                data=mynewspdf2@data)
summary(lm6_res_vi)

#Step_M6_9: Saving coefficients and SE
coef_lm6_res_vi<-summary(margins(lm6_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1700)))$AME
se_lm6_res_vi<-summary(margins(lm6_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1700)))$SE

d6_vi <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d6_vi) <- c("lngdpcap", "p_ottoman", "vienna_distance")
f6_vi <- as.formula("lngdpcap ~ 0 + p_ottoman + vienna_distance")
p6_vi <- lm(f6_vi, d6_vi)
p6_vi$coefficients<-coef_lm6_res_vi
p6_se_vi<- setNames(se_lm6_res_vi, c("p_ottoman", "vienna_distance"))
p6_vi$df.residual<-NA

#Step_M6_10: Run regression with Trade Route Density
lm6_res_tr<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  factor(ctrnum) +
                  route_density_2,
                data=mynewspdf2@data)
summary(lm6_res_tr)

#Step_M6_11: Saving coefficients and SE
coef_lm6_res_tr<-summary(margins(lm6_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1700)))$AME
se_lm6_res_tr<-summary(margins(lm6_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1700)))$SE

d6_tr <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d6_tr) <- c("lngdpcap", "p_ottoman", "route_density_2")
f6_tr <- as.formula("lngdpcap~ 0 + p_ottoman + route_density_2")
p6_tr <- lm(f6_tr, d6_tr)
p6_tr$coefficients<-coef_lm6_res_tr
p6_se_tr<- setNames(se_lm6_res_tr, c("p_ottoman", "route_density_2"))
p6_tr$df.residual<-NA

#Step_M6_12: Spatial Filtering
lm6_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + yearpos+ 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y +
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)


#Step_M6_13: Creating the strings of eigenvalues
mynewspdf3<-mynewspdf2@data
mynewspdf4<-cbind(mynewspdf3, fitted(lm6_res_sp))
eigen_vars<-c(names(mynewspdf4[,grepl( "vec" , names(mynewspdf4))]))

#Step_M6_14: Include fitted results
lm6_res_fitted<-lm(as.formula(
  paste("lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum)+", paste(eigen_vars, collapse = "+"))),
  data=mynewspdf4)
summary(lm6_res_fitted)

#Step_M6_15: Saving coefficients and SE
coef_lm6_res_fitted<-summary(margins(lm6_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$AME
se_lm6_res_fitted<-summary(margins(lm6_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1700)))$SE

d6_fitted <- as.data.frame(matrix(rnorm(nrow(mynewspdf4) * 2), nc = 2))
names(d6_fitted) <- c("lngdpcap", "p_ottoman")
f6_fitted <- as.formula("lngdpcap ~ 0 + p_ottoman")
p6_fitted <- lm(f6_fitted, d6_fitted)
p6_fitted$coefficients<-coef_lm6_res_fitted
p6_fitted_se<- setNames(se_lm6_res_fitted, c("p_ottoman"))
p6_fitted$df.residual<-NA

#Step_M6_16: Save Moran I stats
lm6_moran<-moran.test(x=resid(lm6_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm6_moran_p<-moran.test(x=resid(lm6_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm6_moran_star<-paste(round(lm6_moran, digits = 3), "",
                       ifelse(lm6_moran_p<.01,"***",
                              ifelse(lm6_moran_p<.05,"**",
                                     ifelse(lm6_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 7#
#########

#Step_M7_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                              !is.na(nuts2@data$p_ottoman) &
                              !is.na(nuts2@data$yearpos) &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]


#Step_M7_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M7_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M7_4: Run regression
lm7_res<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum),
             data=mynewspdf2@data)
summary(lm7_res)

#Step_M7_5: Saving coefficients and SE
coef_lm7_res<-summary(margins(lm7_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$AME
se_lm7_res<-summary(margins(lm7_res,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$SE

d7 <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 2), nc = 2))
names(d7) <- c("lngdpcap", "p_ottoman")
f7 <- as.formula("lngdpcap ~ 0 + p_ottoman")
p7 <- lm(f7, d7)
p7$coefficients<-coef_lm7_res
p7_se<- setNames(se_lm7_res, c("p_ottoman"))
p7$df.residual<-NA

#Step_M7_6: Run regression with distance to Istanbul
lm7_res_di<-lm(lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               istanbul_distance+
               factor(ctrnum),
             data=mynewspdf2@data)
summary(lm7_res_di)

#Step_M7_7: Saving coefficients and SE
coef_lm7_res_di<-summary(margins(lm7_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1750)))$AME
se_lm7_res_di<-summary(margins(lm7_res_di,  vce = "delta", variables = c("p_ottoman", "istanbul_distance"), at = list(yearpos = 1750)))$SE

d7_di <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d7_di) <- c("lngdpcap", "p_ottoman", "istanbul_distance")
f7_di <- as.formula("lngdpcap~ 0 + p_ottoman + istanbul_distance")
p7_di <- lm(f7_di, d7_di)
p7_di$coefficients<-coef_lm7_res_di
p7_se_di<- setNames(se_lm7_res_di, c("p_ottoman", "istanbul_distance"))
p7_di$df.residual<-NA

#Step_M7_8: Run regression with distance to Vienna
lm7_res_vi<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  vienna_distance+
                 factor(ctrnum),
                data=mynewspdf2@data)
summary(lm7_res_vi)

#Step_M7_9: Saving coefficients and SE
coef_lm7_res_vi<-summary(margins(lm7_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1750)))$AME
se_lm7_res_vi<-summary(margins(lm7_res_vi,  vce = "delta", variables = c("p_ottoman", "vienna_distance"), at = list(yearpos = 1750)))$SE

d7_vi <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d7_vi) <- c("lngdpcap", "p_ottoman", "vienna_distance")
f7_vi <- as.formula("lngdpcap ~ 0 + p_ottoman + vienna_distance")
p7_vi <- lm(f7_vi, d7_vi)
p7_vi$coefficients<-coef_lm7_res_vi
p7_se_vi<- setNames(se_lm7_res_vi, c("p_ottoman", "vienna_distance"))
p7_vi$df.residual<-NA

#Step_M7_10: Run regression with Trade Route Density
lm7_res_tr<-lm(lngdpcap~p_ottoman * yearpos+ 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  route_density_2 +
                 factor(ctrnum),
                data=mynewspdf2@data)
summary(lm7_res_tr)

#Step_M7_11: Saving coefficients and SE
coef_lm7_res_tr<-summary(margins(lm7_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1750)))$AME
se_lm7_res_tr<-summary(margins(lm7_res_tr,  vce = "delta", variables = c("p_ottoman", "route_density_2"), at = list(yearpos = 1750)))$SE

d7_tr <- as.data.frame(matrix(rnorm(nrow(mynewspdf2@data) * 3), nc = 3))
names(d7_tr) <- c("lngdpcap", "p_ottoman", "route_density_2")
f7_tr <- as.formula("lngdpcap ~ 0 + p_ottoman + route_density_2")
p7_tr <- lm(f7_tr, d7_tr)
p7_tr$coefficients<-coef_lm7_res_tr
p7_se_tr<- setNames(se_lm7_res_tr, c("p_ottoman", "route_density_2"))
p7_tr$df.residual<-NA

#Step_M7_12: Spatial Filtering
lm7_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + yearpos+ 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y +
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M7_13: Creating the strings of eigenvalues
mynewspdf3<-mynewspdf2@data
mynewspdf4<-cbind(mynewspdf3, fitted(lm7_res_sp))
eigen_vars<-c(names(mynewspdf4[,grepl( "vec" , names(mynewspdf4))]))

#Step_M7_14: Include fitted results
lm7_res_fitted<-lm(as.formula(
  paste("lngdpcap~p_ottoman * yearpos+ 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y +
               factor(ctrnum)+", paste(eigen_vars, collapse = "+"))),
  data=mynewspdf4)
summary(lm7_res_fitted)

#Step_M7_15: Saving coefficients and SE
coef_lm7_res_fitted<-summary(margins(lm7_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$AME
se_lm7_res_fitted<-summary(margins(lm7_res_fitted,  vce = "delta", variables = "p_ottoman", at = list(yearpos = 1750)))$SE

d7_fittted <- as.data.frame(matrix(rnorm(nrow(mynewspdf4) * 2), nc = 2))
names(d7_fittted) <- c("lngdpcap", "p_ottoman")
f7_fitted <- as.formula("lngdpcap ~ 0 + p_ottoman")
p7_fitted <- lm(f7_fitted, d7_fittted)
p7_fitted$coefficients<-coef_lm7_res_fitted
p7_fitted_se<- setNames(se_lm7_res_fitted, c("p_ottoman"))
p7_fitted$df.residual<-NA

#Step_M7_16: Save Moran I stats
lm7_moran<-moran.test(x=resid(lm7_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm7_moran_p<-moran.test(x=resid(lm7_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm7_moran_star<-paste(round(lm7_moran, digits = 3), "",
                       ifelse(lm7_moran_p<.01,"***",
                              ifelse(lm7_moran_p<.05,"**",
                                     ifelse(lm7_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 8#
#########

#Step_M8_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$firmscap) &
                              !is.na(nuts2@data$p_ottoman)  &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]

#Step_M8_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M8_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M8_4: Run regression
felm8_res<-felm(firmscap~p_ottoman + 
              tavg + lograin + logelev +
              logslope + wheat_rainfed + log_sea_lines_distance +
              log_large_rivers_distance + point_x + point_y 
              |ctrnum | 0 | ctrnum,
            data=mynewspdf2@data)
summary(felm8_res)
#Saving Cluster SE
felm8_res_cse<-felm8_res$cse

#Step_M8_5: Run regression with distance to Istanbul
felm8_res_di<-felm(firmscap~p_ottoman + 
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y +
                 istanbul_distance 
                 |ctrnum | 0 |ctrnum,
               data=mynewspdf2@data)
summary(felm8_res_di)
#Saving Cluster SE
felm8_res_cse_di<-felm8_res_di$cse

#Step_M8_6: Run regression with distance to Vienna
felm8_res_vi<-felm(firmscap~p_ottoman + 
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y +
                 vienna_distance 
                 |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm8_res_vi)
#Saving Cluster SE
felm8_res_cse_vi<-felm8_res_vi$cse

#Step_M8_7: Run regression with Trade Route Density
felm8_res_tr<-felm(firmscap~p_ottoman + route_density_2+
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y
                 |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm8_res_tr)
#Saving Cluster SE
felm8_res_cse_tr<-felm8_res_tr$cse

#Step_M8_8: Spatial Filtering
lm8_res_sp<-SpatialFiltering(firmscap~p_ottoman + 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y+
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M8_9: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm8_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M8_10: Include fitted results
x<-as.formula(paste("firmscap~p_ottoman + 
                      tavg + lograin + logelev +
                        logslope + wheat_rainfed + log_sea_lines_distance +
                        log_large_rivers_distance + point_x + point_y +", paste(eigen_string), 
                      "|ctrnum | 0 | ctrnum"))

felm8_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm8_res_fitted)
#Saving Cluster SE
felm8_res_fitted_cse<-felm8_res_fitted$cse

#Step_M8_11: Save Moran I stats
lm8_moran<-moran.test(x=resid(felm8_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm8_moran_p<-moran.test(x=resid(felm8_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm8_moran_star<-paste(round(lm8_moran, digits = 3), "",
                       ifelse(lm8_moran_p<.01,"***",
                              ifelse(lm8_moran_p<.05,"**",
                                     ifelse(lm8_moran_p<.1,"*",
                                            ""))), "", "", sep = "")


#########
#MODEL 9#
#########

#Step_M9_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$firms10cap) &
                              !is.na(nuts2@data$p_ottoman)  &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]

#Step_M9_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M9_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M9_4: Run regression
felm9_res<-felm(firms10cap~p_ottoman + 
              tavg + lograin + logelev +
              logslope + wheat_rainfed + log_sea_lines_distance +
              log_large_rivers_distance + point_x + point_y
              |ctrnum | 0 | ctrnum,
            data=mynewspdf2@data)
summary(felm9_res)
#Saving Cluster SE
felm9_res_cse<-felm9_res$cse

#Step_M9_5: Run regression with distance to Istanbul
felm9_res_di<-felm(firms10cap~p_ottoman + 
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y +
                 istanbul_distance 
                 |ctrnum | 0 |ctrnum,
               data=mynewspdf2@data)
summary(felm9_res_di)
#Saving Cluster SE
felm9_res_cse_di<-felm9_res_di$cse

#Step_M9_6: Run regression with distance to Vienna
felm9_res_vi<-felm(firms10cap~p_ottoman + 
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y +
                 vienna_distance 
                 |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm9_res_vi)
#Saving Cluster SE
felm9_res_cse_vi<-felm9_res_vi$cse

#Step_M9_7: Run regression with Trade Route Density
felm9_res_tr<-felm(firms10cap~p_ottoman + route_density_2+
                 tavg + lograin + logelev +
                 logslope + wheat_rainfed + log_sea_lines_distance +
                 log_large_rivers_distance + point_x + point_y 
                 |ctrnum | 0 | ctrnum,
               data=mynewspdf2@data)
summary(felm9_res_tr)
#Saving Cluster SE
felm9_res_cse_tr<-felm9_res_tr$cse

#Step_M9_8: Spatial Filtering
lm9_res_sp<-SpatialFiltering(firms10cap~p_ottoman + 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y+
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M9_9: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm9_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M9_10: Include fitted results
x<-as.formula(paste("firms10cap~p_ottoman + 
                        tavg + lograin + logelev +
                        logslope + wheat_rainfed + log_sea_lines_distance +
                        log_large_rivers_distance + point_x + point_y + ", paste(eigen_string),
                      "|ctrnum | 0 | ctrnum"))
felm9_res_fitted<-felm(x,
                   data=mynewspdf3)
summary(felm9_res_fitted)
#Saving Cluster SE
felm9_res_fitted_cse<-felm9_res_fitted$cse

#Step_M9_11: Save Moran I stats
lm9_moran<-moran.test(x=resid(felm9_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm9_moran_p<-moran.test(x=resid(felm9_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm9_moran_star<-paste(round(lm9_moran, digits = 3), "",
                       ifelse(lm9_moran_p<.01,"***",
                              ifelse(lm9_moran_p<.05,"**",
                                     ifelse(lm9_moran_p<.1,"*",
                                            ""))), "", "", sep = "")

##########
#MODEL 10#
##########

#Step_M10_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lntmapplic) &
                              !is.na(nuts2@data$p_ottoman)  &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]


#Step_M10_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M10_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M10_4: Run regression
felm10_res<-felm(lntmapplic~p_ottoman + 
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y
               |ctrnum | 0 | ctrnum,
             data=mynewspdf2@data)
summary(felm10_res)
#Saving Cluster SE
felm10_res_cse<-felm10_res$cse

#Step_M10_5: Run regression with distance to Istanbul
felm10_res_di<-felm(lntmapplic~p_ottoman + 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  istanbul_distance 
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm10_res_di)
#Saving Cluster SE
felm10_res_cse_di<-felm10_res_di$cse

#Step_M10_6: Run regression with distance to Vienna
felm10_res_vi<-felm(lntmapplic~p_ottoman + 
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  vienna_distance 
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm10_res_vi)
#Saving Cluster SE
felm10_res_cse_vi<-felm10_res_vi$cse

#Step_M10_7: Run regression with Trade Route Density
felm10_res_tr<-felm(lntmapplic~p_ottoman + route_density_2+
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y 
                |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm10_res_tr)
#Saving Cluster SE
felm10_res_cse_tr<-felm10_res_tr$cse

#Step_M10_8: Spatial Filtering
lm10_res_sp<-SpatialFiltering(lntmapplic~p_ottoman + 
                                tavg + lograin + logelev +
                                logslope + wheat_rainfed + log_sea_lines_distance +
                                log_large_rivers_distance + point_x + point_y+
                                factor(ctrnum),
                              data = mynewspdf2, 
                              nb = list_nb, 
                              style="W", 
                              ExactEV = F, zero.policy = T)

#Step_M10_9: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm10_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M10_10: Include fitted results
x<-as.formula(paste("lntmapplic~p_ottoman + 
                        tavg + lograin + logelev +
                        logslope + wheat_rainfed + log_sea_lines_distance +
                        log_large_rivers_distance + point_x + point_y + ", paste(eigen_string),
                      "|ctrnum | 0 | ctrnum"))
felm10_res_fitted<-felm(x,
                    data=mynewspdf3)
summary(felm10_res_fitted)
#Saving Cluster SE
felm10_res_fitted_cse<-felm10_res_fitted$cse

#Step_M10_11: Save Moran I stats
lm10_moran<-moran.test(x=resid(felm10_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm10_moran_p<-moran.test(x=resid(felm10_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm10_moran_star<-paste(round(lm10_moran, digits = 3), "",
                       ifelse(lm10_moran_p<.01,"***",
                              ifelse(lm10_moran_p<.05,"**",
                                     ifelse(lm10_moran_p<.1,"*",
                                            ""))), "", "", sep = "")

###########
#MODEL 11#
##########

#Step_M11_1: Remove NAs
mynewspdf2<-nuts2[which(!is.na(nuts2@data$lngdpcap) &
                              !is.na(nuts2@data$lntmapplic) &
                              !is.na(nuts2@data$p_ottoman)  &
                              !is.na(nuts2@data$ctrnum) &
                              !is.na(nuts2@data$tavg) &
                              !is.na(nuts2@data$lograin) &
                              !is.na(nuts2@data$logelev) &
                              !is.na(nuts2@data$logslope) & 
                              !is.na(nuts2@data$wheat_rainfed) &
                              !is.na(nuts2@data$log_sea_lines_distance) & 
                              !is.na(nuts2@data$log_large_rivers_distance) & 
                              !is.na(nuts2@data$point_x) & 
                              !is.na(nuts2@data$point_y)),]


#Step_M11_2: Create list of neighbors
list_nb=poly2nb(mynewspdf2, row.names = mynewspdf2$merge_id)

#Step_M11_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M11_4: Run regression
felm11_res<-felm(lngdpcap~p_ottoman + lntmapplic+
               tavg + lograin + logelev +
               logslope + wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + point_x + point_y 
               |ctrnum | 0 | ctrnum,
             data=mynewspdf2@data)
summary(felm11_res)
#Saving Cluster SE
felm11_res_cse<-felm11_res$cse

#Step_M11_5: Run regression with distance to Istanbul
felm11_res_di<-felm(lngdpcap~p_ottoman + lntmapplic+
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  istanbul_distance
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm11_res_di)
#Saving Cluster SE
felm11_res_cse_di<-felm11_res_di$cse

#Step_M11_6: Run regression with distance to Vienna
felm11_res_vi<-felm(lngdpcap~p_ottoman + lntmapplic+
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y +
                  vienna_distance
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm11_res_vi)
#Saving Cluster SE
felm11_res_cse_vi<-felm11_res_vi$cse

#Step_M11_7: Run regression with Trade Route Density
felm11_res_tr<-felm(lngdpcap~p_ottoman + lntmapplic+ route_density_2+
                  tavg + lograin + logelev +
                  logslope + wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + point_x + point_y 
                |ctrnum | 0 | ctrnum,
                data=mynewspdf2@data)
summary(felm11_res_tr)
#Saving Cluster SE
felm11_res_cse_tr<-felm11_res_tr$cse

#Step_M11_8: Spatial Filtering
lm11_res_sp<-SpatialFiltering(lngdpcap~p_ottoman + lntmapplic+
                                 tavg + lograin + logelev +
                                 logslope + wheat_rainfed + log_sea_lines_distance +
                                 log_large_rivers_distance + point_x + point_y +
                                 factor(ctrnum),
                               data = mynewspdf2, 
                               nb = list_nb, 
                               style="W", 
                               ExactEV = F, zero.policy = T)

#Step_M11_9: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm11_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf2@data, new_df)

#Step_M11_10: Include fitted results
x<-as.formula(paste("lngdpcap~p_ottoman + lntmapplic+
                         tavg + lograin + logelev +
                         logslope + wheat_rainfed + log_sea_lines_distance +
                         log_large_rivers_distance + point_x + point_y + ", paste(eigen_string),
                      "|ctrnum | 0 | ctrnum"))
felm11_res_fitted<-felm(x,
                    data=mynewspdf3)
summary(felm11_res_fitted)
#Saving Cluster SE
felm11_res_fitted_cse<-felm11_res_fitted$cse

#Step_M11_11: Save Moran I stats
lm11_moran<-moran.test(x=resid(felm11_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm11_moran_p<-moran.test(x=resid(felm11_res_fitted),listw=list_nb_w, zero.policy = T)$p.value

lm11_moran_star<-paste(round(lm11_moran, digits = 3), "",
                        ifelse(lm11_moran_p<.01,"***",
                               ifelse(lm11_moran_p<.05,"**",
                                      ifelse(lm11_moran_p<.1,"*",
                                             ""))), "", "", sep = "")


################
#MODEL 12 FILES#
################

#Step1: Reading DTA files
educ_df<-read_dta("./data/educ_data.dta")

#Step2: Reading shapefile for NUTS polygons (counties)
nuts_sh <- readOGR(dsn="./data/shape_files.gdb",
                   layer="nuts_EU_2001", use_iconv = TRUE, encoding = "UTF-8")
#Step3: Merge
nuts_edu<-merge(nuts_sh, educ_df, by.x = "MERGE_ID_EDU", by.y = "MERGE_ID_EDU")


##########
#MODEL 12#
##########
#Step_M12_1: Remove NAs
mynewspdf_2<-nuts_edu[which(!is.na(nuts_edu@data$secondaryed2) &
                                 !is.na(nuts_edu@data$p_ottoman) &
                                 !is.na(nuts_edu@data$route_density_2) &
                                 !is.na(nuts_edu@data$ctrnum) &
                                 !is.na(nuts_edu@data$tavg) &
                                 !is.na(nuts_edu@data$lograin) &
                                 !is.na(nuts_edu@data$logelev) &
                                 !is.na(nuts_edu@data$logslope) & 
                                 !is.na(nuts_edu@data$point_x) & 
                                 !is.na(nuts_edu@data$point_y) &
                                 !is.na(nuts_edu@data$wheat_rainfed) &
                                 !is.na(nuts_edu@data$log_sea_lines_distance) &
                                 !is.na(nuts_edu@data$log_large_rivers_distance)),]


#Step_M12_2: Create list of neighbors
list_nb=poly2nb(mynewspdf_2, row.names = mynewspdf_2$MERGE_ID_EDU)

#Step_M12_3: Convert nb to listw type
list_nb_w=nb2listw(list_nb, zero.policy=TRUE)
print(list_nb_w, zero.policy=TRUE)

#Step_M12_4: Run regression
felm12_res<-felm(secondaryed2~p_ottoman +
               tavg + lograin + logelev +
               logslope + 
               wheat_rainfed + log_sea_lines_distance +
               log_large_rivers_distance + 
                 point_x + point_y 
               |ctrnum | 0 | ctrnum,
             data=mynewspdf_2@data)
summary(felm12_res)
#Saving Cluster SE
felm12_res_cse<-felm12_res$cse

#Step_M12_5: Run regression with distance to Istanbul
felm12_res_di<-felm(secondaryed2~p_ottoman +
                  tavg + lograin + logelev +
                  logslope + 
                  wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + 
                  #  point_x + point_y +
                  istanbul_distance 
                  |ctrnum | 0 |ctrnum,
                data=mynewspdf_2@data)
summary(felm12_res_di)
#Saving Cluster SE
felm12_res_cse_di<-felm12_res_di$cse

#Step_M12_6: Run regression with distance to Vienna
felm12_res_vi<-felm(secondaryed2~p_ottoman +
                  tavg + lograin + logelev +
                  logslope + 
                  wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + 
                  #point_x + point_y
                  vienna_distance 
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf_2@data)
summary(felm12_res_vi)
#Saving Cluster SE
felm12_res_cse_vi<-felm12_res_vi$cse

#Step_M12_7: Run regression with Trade Route Density
felm12_res_tr<-felm(secondaryed2~p_ottoman + route_density_2+
                  tavg + lograin + logelev +
                  logslope + 
                  wheat_rainfed + log_sea_lines_distance +
                  log_large_rivers_distance + 
                    point_x + point_y
                  |ctrnum | 0 | ctrnum,
                data=mynewspdf_2@data)
summary(felm12_res_tr)
#Saving Cluster SE
felm12_res_cse_tr<-felm12_res_tr$cse

#Step_M12_8: Spatial Filtering
lm12_res_sp<-SpatialFiltering(secondaryed2~p_ottoman +
                                 tavg + lograin + logelev +
                                 logslope + 
                                 wheat_rainfed + log_sea_lines_distance +
                                 log_large_rivers_distance + 
                                 #Including Lat and Lon causes multicollinearity
                                 #point_x + point_y
                                 factor(ctrnum),
                               data = mynewspdf_2, 
                               nb = list_nb, 
                               style="W", 
                               ExactEV = F, zero.policy = T)

#Step_M12_9: Creating the strings of eigenvalues
new_df<-as.data.frame(fitted(lm12_res_sp))
eigen_vec<-names(new_df)
eigen_string =  paste(eigen_vec, collapse=' + ' )
mynewspdf3<-cbind(mynewspdf_2@data, new_df)

#Step_M12_10: Include fitted results
x<-as.formula(paste("secondaryed2~p_ottoman +
                         istanbul_distance+
                         tavg + lograin + logelev +
                         logslope +
                         point_x + point_y + ", paste(eigen_string), 
                      "|ctrnum | 0 | ctrnum"))
felm12_res_fitted<-felm(x,
                    data=mynewspdf3)
summary(felm12_res_fitted)
#Saving Cluster SE
felm12_res_fitted_cse<-felm12_res_fitted$cse

#Step_M12_11: Save Moran I stats
lm12_moran<-moran.test(x=resid(felm12_res_fitted),listw=list_nb_w, zero.policy = T)$statistic
lm12_moran_p<-moran.test(x=resid(felm12_res_fitted),listw=list_nb_w, zero.policy = T)$p.value
lm12_moran_star<-paste(round(lm12_moran, digits = 3), "",
                        ifelse(lm12_moran_p<.01,"***",
                               ifelse(lm12_moran_p<.05,"**",
                                      ifelse(lm12_moran_p<.1,"*",
                                             ""))), "", "", sep = "")
                               


####################################################
#MAIN TEXT Results with no Sp. Autocorr - TABLE1 FELM#
####################################################

table1 = stargazer(felm1_res,
                           felm2_res,
                           felm3_res,
                           felm4_res,
                           p5,
                           p6,
                           p7,
                           felm8_res,
                           felm9_res,
                           felm10_res,
                           felm11_res,
                           felm12_res,
                           se = list(felm1_res_cse,
                                     felm2_res_cse,
                                     felm3_res_cse,
                                     felm4_res_cse,
                                     p5_se,
                                     p6_se,
                                     p7_se,
                                     felm8_res_cse,
                                     felm9_res_cse,
                                     felm10_res_cse,
                                     felm11_res_cse,
                                     felm12_res_cse),
                           column.labels= c(rep("\\shortstack{Log\\\\GDP/cap}", 7),
                                            "\\shortstack{Firms\\\\per\\\\Capita}",
                                            "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
                                            "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
                                            "\\shortstack{Log\\\\GDP/cap}",
                                            "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}"),
                           keep = c("p_ottoman",
                                    "route_density_2", 
                                    "lntmapplic"),
                           covariate.labels=c("Prop. Ottoman",
                                              "Trade Route Dens.", 
                                              "Log Trade Mark Certif. per Capita"),
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           dep.var.caption = "",
                           dep.var.labels.include = F,
                           model.names=F,
                           omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                           add.lines=list(c("Country FE", c("No", rep("Yes", 12))),
                                          c("Model", rep("OLS", 12)),
                                          c("Controls", c("No", "No", rep("Geo + Pos", 10))),
                                          c("Years", c(rep("All", 4), "1650", "1700", "1750", rep("All", 5)))))
note_latex <- "\\multicolumn{13}{l} {\\parbox[t]{38cm}{ \\footnotesize \\textit{Notes:} 
All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level.
All models are estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude. 
Models 5, 6, and 7 present results from models in which
proportion Ottoman is interacted with the average year of Ottoman rule, 
and marginal effects in 1650, 1700, and 1750 are calculated.
*** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
table1 [grepl("Note", table1)] <- note_latex
cat(table1, file="./Paper/tables/table1.tex", sep="\n")


###############################
#CLUSTER RESULTS FELM + SP. AUTO#
###############################
tableA3 = stargazer(felm1_res_fitted,
                          felm2_res_fitted,
                          felm3_res_fitted,
                          felm4_res_fitted,
                          p5_fitted,
                          p6_fitted,
                          p7_fitted,
                          felm8_res_fitted,
                          felm9_res_fitted,
                          felm10_res_fitted,
                          felm11_res_fitted,
                          felm12_res_fitted,
                          se = list(felm1_res_fitted_cse,
                                    felm2_res_fitted_cse,
                                    felm3_res_fitted_cse,
                                    felm4_res_fitted_cse,
                                    p5_fitted_se,
                                    p6_fitted_se,
                                    p7_fitted_se,
                                    felm8_res_fitted_cse,
                                    felm9_res_fitted_cse,
                                    felm10_res_fitted_cse,
                                    felm11_res_fitted_cse,
                                    felm12_res_fitted_cse),
                          column.labels= c(rep("\\shortstack{Log\\\\GDP/cap}", 7),
                                           "\\shortstack{Firms\\\\per\\\\Capita}",
                                           "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
                                           "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
                                           "\\shortstack{Log\\\\GDP/cap}",
                                           "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}"),
                          keep = c("p_ottoman",
                                   "route_density_2", 
                                   "lntmapplic"),
                          covariate.labels=c("Prop. Ottoman",
                                             "Trade Route Dens.", 
                                             "Log Trade Mark Certif. per Capita"),
                          star.char = c("+", "*", "**", "***"),
                          star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                          dep.var.caption = "",
                          dep.var.labels.include = F,
                          model.names=F,
                          omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                          add.lines=list(c("Country FE", c("No", rep("Yes", 12))),
                                         c("Model", rep("OLS", 12)),
                                         c("Controls", c("No", "No", rep("Geo + Pos", 10))),
                                         c("Years", c(rep("All", 4), "1650", "1700", "1750", rep("All", 5))),
                                         c("Moran eigenvectors", c(rep("Yes", 12))),
                                         c("Moran'I residual", c(lm1_moran_star,
                                                                 lm2_moran_star,
                                                                 lm3_moran_star,
                                                                 lm4_moran_star,
                                                                 lm5_moran_star,
                                                                 lm6_moran_star,
                                                                 lm7_moran_star,
                                                                 lm8_moran_star,
                                                                 lm9_moran_star,
                                                                 lm10_moran_star,
                                                                 lm11_moran_star,
                                                                 lm12_moran_star))))
note_latex <- "\\multicolumn{13}{l} {\\parbox[t]{38cm}{ \\footnotesize \\textit{Notes:} 
All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level.
All models are estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude. 
Models 5, 6, and 7 present results from models in which
proportion Ottoman is interacted with the average year of Ottoman rule, 
and marginal effects in 1650, 1700, and 1750 are calculated.
*** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
tableA3 [grepl("Note", tableA3)] <- note_latex
#cat(pct_no_water, file="pct_no_water.tex", sep="\n")
cat(tableA3, file="./Paper/tables/tableA3.tex", sep="\n")


############################
#Distance to Istanbul FELM#
###########################

tableA8 = stargazer(felm1_res_di,
                           felm2_res_di,
                           felm3_res_di,
                           felm4_res_di,
                           p5_di,
                           p6_di,
                           p7_di,
                           felm8_res_di,
                           felm9_res_di,
                           felm10_res_di,
                           felm11_res_di,
                           felm12_res_di,
                           se = list(felm1_res_cse_di,
                                     felm2_res_cse_di,
                                     felm3_res_cse_di,
                                     felm4_res_cse_di,
                                     p5_se_di,
                                     p6_se_di,
                                     p7_se_di,
                                     felm8_res_cse_di,
                                     felm9_res_cse_di,
                                     felm10_res_cse_di,
                                     felm11_res_cse_di,
                                     felm12_res_cse_di),
                           column.labels= c(rep("\\shortstack{Log\\\\GDP/cap}", 7),
                                            "\\shortstack{Firms\\\\per\\\\Capita}",
                                            "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
                                            "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
                                            "\\shortstack{Log\\\\GDP/cap}",
                                            "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}"),
                           keep = c("p_ottoman", 
                                    "route_density_2", 
                                    "lntmapplic",
                                    "istanbul_distance"),
                           covariate.labels=c("Prop. Ottoman", 
                                              "Trade Route Dens.", 
                                              "Log Trade Mark Certif. per Capita",
                                              "Istanbul Distance"),
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           dep.var.caption = "",
                           dep.var.labels.include = F,
                           model.names=F,
                           omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                           add.lines=list(c("Country FE", c("No", rep("Yes", 12))),
                                          c("Model", rep("OLS", 12)),
                                          c("Controls", c("No", "No", rep("Geo + Pos", 10))),
                                          c("Years", c(rep("All", 4), "1650", "1700", "1750", rep("All", 5)))))
note_latex <- "\\multicolumn{13}{l} {\\parbox[t]{38cm}{ \\footnotesize \\textit{Notes:} 
All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level.
All models are estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude. 
Models 5, 6, and 7 present results from models in which
proportion Ottoman is interacted with the average year of Ottoman rule, 
and marginal effects in 1650, 1700, and 1750 are calculated.
*** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
tableA8[grepl("Note", tableA8)] <- note_latex
cat(tableA8, file="./Paper/tables/tableA8.tex", sep="\n")


#########################
#Distance to VIENNA FELM#
#########################

tableA9 = stargazer(felm1_res_vi,
                           felm2_res_vi,
                           felm3_res_vi,
                           felm4_res_vi,
                           p5_vi,
                           p6_vi,
                           p7_vi,
                           felm8_res_vi,
                           felm9_res_vi,
                           felm10_res_vi,
                           felm11_res_vi,
                           felm12_res_vi,
                           se = list(felm1_res_cse_vi,
                                     felm2_res_cse_vi,
                                     felm3_res_cse_vi,
                                     felm4_res_cse_vi,
                                     p5_se_vi,
                                     p6_se_vi,
                                     p7_se_vi,
                                     felm8_res_cse_vi,
                                     felm9_res_cse_vi,
                                     felm10_res_cse_vi,
                                     felm11_res_cse_vi,
                                     felm12_res_cse_vi),
                           column.labels= c(rep("\\shortstack{Log\\\\GDP/cap}", 7),
                                            "\\shortstack{Firms\\\\per\\\\Capita}",
                                            "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
                                            "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
                                            "\\shortstack{Log\\\\GDP/cap}",
                                            "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}"),
                           keep = c("p_ottoman", 
                                    "route_density_2", 
                                    "lntmapplic",
                                    "vienna_distance"),
                           covariate.labels=c("Prop. Ottoman", 
                                              "Trade Route Dens.", 
                                              "Log Trade Mark Certif. per Capita",
                                              "Vienna Distance"),
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           dep.var.caption = "",
                           dep.var.labels.include = F,
                           model.names=F,
                           omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                           add.lines=list(c("Country FE", c("No", rep("Yes", 12))),
                                          c("Model", rep("OLS", 12)),
                                          c("Controls", c("No", "No", rep("Geo + Pos", 10))),
                                          c("Years", c(rep("All", 4), "1650", "1700", "1750", rep("All", 5)))))
note_latex <- "\\multicolumn{13}{l} {\\parbox[t]{38cm}{ \\footnotesize \\textit{Notes:} 
All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level.
All models are estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude. 
Models 5, 6, and 7 present results from models in which
proportion Ottoman is interacted with the average year of Ottoman rule, 
and marginal effects in 1650, 1700, and 1750 are calculated.
*** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
tableA9[grepl("Note", tableA9)] <- note_latex
cat(tableA9, file="./Paper/tables/tableA9.tex", sep="\n")


##########################
#CONTROL FOR TRADE ROUTES#
##########################

tableA10 = stargazer(felm1_res,
                           felm2_res,
                           felm3_res,
                           felm4_res,
                           p5_tr,
                           p6_tr,
                           p7_tr,
                           felm8_res_tr,
                           felm9_res_tr,
                           felm10_res_tr,
                           felm11_res_tr,
                           felm12_res_tr,
                           se = list(felm1_res_cse,
                                     felm2_res_cse,
                                     felm3_res_cse,
                                     felm4_res_cse,
                                     p5_se_tr,
                                     p6_se_tr,
                                     p7_se_tr,
                                     felm8_res_cse_tr,
                                     felm9_res_cse_tr,
                                     felm10_res_cse_tr,
                                     felm11_res_cse_tr,
                                     felm12_res_cse_tr),
                           column.labels= c(rep("\\shortstack{Log\\\\GDP/cap}", 7),
                                            "\\shortstack{Firms\\\\per\\\\Capita}",
                                            "\\shortstack{Firms\\\\> 10 employees\\\\per Capita}",
                                            "\\shortstack{Trade Mark\\\\Certifications\\\\per Capita}",
                                            "\\shortstack{Log\\\\GDP/cap}",
                                            "\\shortstack{Share. Pop.\\\\with Sec.\\\\Education}"),
                           keep = c("p_ottoman",
                                    "lntmapplic",
                                    "route_density_2"),
                           covariate.labels=c("Prop. Ottoman",
                                              "Log Trade Mark Certif. per Capita",
                                              "Trade Route Dens."),
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           dep.var.caption = "",
                           dep.var.labels.include = F,
                           model.names=F,
                           omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                           add.lines=list(c("Country FE", c("No", rep("Yes", 12))),
                                          c("Model", rep("OLS", 12)),
                                          c("Controls", c("No", "No", rep("Geo + Pos", 10))),
                                          c("Years", c(rep("All", 4), "1650", "1700", "1750", rep("All", 5)))))
note_latex <- "\\multicolumn{13}{l} {\\parbox[t]{38cm}{ \\footnotesize \\textit{Notes:} 
All columns present the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 or on the other dependent variables at the regional level.
All models are estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude. 
Models 5, 6, and 7 present results from models in which
proportion Ottoman is interacted with the average year of Ottoman rule, 
and marginal effects in 1650, 1700, and 1750 are calculated.
*** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
tableA10 [grepl("Note", tableA10)] <- note_latex
cat(tableA10, file="./Paper/tables/tableA10.tex", sep="\n")


#####################
#Effect of Communism#
#####################

df_comm<-mynewspdf2[which(!is.na(mynewspdf2@data$lngdpcap) &
                              !is.na(mynewspdf2@data$p_ottoman) &
                              !is.na(mynewspdf2@data$ctrnum) &
                              !is.na(mynewspdf2@data$tavg) &
                              !is.na(mynewspdf2@data$lograin) &
                              !is.na(mynewspdf2@data$logelev) &
                              !is.na(mynewspdf2@data$logslope) & 
                              !is.na(mynewspdf2@data$wheat_rainfed) &
                              !is.na(mynewspdf2@data$log_sea_lines_distance) & 
                              !is.na(mynewspdf2@data$log_large_rivers_distance) & 
                              !is.na(mynewspdf2@data$point_x) & 
                              !is.na(mynewspdf2@data$point_y)),]


lm3_res_di_com<-lm(lngdpcap~p_ottoman + communist +
                        tavg + lograin + logelev +
                        logslope + wheat_rainfed + log_sea_lines_distance +
                        log_large_rivers_distance + point_x + point_y + factor(ctrnum),
                      data=df_comm@data)
summary(lm3_res_di_com)

#Function for Cluster SE
se_robust_cluster <- function(x){
  res<-coeftest(x, vcov = cluster.vcov(x, df_comm@data$ctrnum, force_posdef=T))[, 2]
  return(res)}

lm3_res_se_di_com<-se_robust_cluster(lm3_res_di_com)
  

tableA5 = stargazer(lm3_res_di_com,
                           se = list(lm3_res_se_di_com),
                           column.labels= c("\\shortstack{Log\\\\GDP/cap}"),
                           keep = c("p_ottoman", 
                                    "communist"),
                           covariate.labels=c("Prop. Ottoman", 
                                              "Communist Past"),
                           star.char = c("+", "*", "**", "***"),
                           star.cutoffs = c(0.1, 0.05, 0.01, 0.001),
                           dep.var.caption = "",
                           dep.var.labels.include = F,
                           model.names=F,
                           omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                           add.lines=list(c("Country FE", c("Yes", "Yes")),
                                          c("Model", rep("OLS", 12)),
                                          c("Controls", c("Geo + Pos")),
                                          c("Years", c(rep("All", 4)))))
note_latex <- "\\multicolumn{2}{l} {\\parbox[t]{17cm}{ \\footnotesize \\textit{Notes:} 
The column presents the marginal effects of the proportion Ottoman indicator (range 0-1) on ln(GDP/capita) in 2016 at the regional level.
The model is estimated with standard errors clustered at the country level. 
Geo controls include average temperature, ln(average altitude), ln(average slope), ln(average precipitation),
wheat growth suitability, distance from rivers, and distance from the sea. 
Pos controls include latitude and longitude.
  *** = p<.001; ** = p<.01, *=p<.05, +=p<.10.}}"
tableA5[grepl("Note", tableA5)] <- note_latex
cat(tableA5, file="./Paper/tables/tableA5.tex", sep="\n")

